## Multiple imputation of welfare spending ##
rm(list=ls())
library(foreign)
library(Amelia)
library(Zelig)

# Load dataset
file <- "Data for imputation.dta"
data <- read.dta(file)
summary(data)

junior.data <- subset(data, select=-c(senior_spell, senior_spell2, senior_spell3, lpolity2, lucdp_interstate, lpurges1))
junior.data <- subset(junior.data, year>=1950)
summary(junior.data)

senior.data <- subset(data, select=-c(junior_spell, junior_spell2, junior_spell3, lpolity2, lucdp_interstate, lpurges1))
senior.data <- subset(senior.data, year>=1950)
summary(senior.data)

# Impute data
junior.imputations <- amelia(junior.data, m = 50, 
                      idvars = NULL, 
                      splinetime = NULL,
                      polytime=3,
                      priors = NULL,
                      intercs = FALSE, 
                      empri = 0.001*nrow(junior.data),
                      ts = "year", 
                      cs = "cowcode", 
                      lags = c("lwelfg", "lln_milex_pc", "lmilex_exp", "lln_gdp", "lln_rents"), 
                      leads = c("lwelfg", "lln_milex_pc", "lmilex_exp", "lln_gdp", "lln_rents"), 
                      logs = NULL, 
                      ords = NULL, 
                      noms = c("juniorcoup_Holger", "seniorcoup_Holger"), 
                      bounds = NULL)

senior.imputations <- amelia(senior.data, m = 50, 
                             idvars = NULL, 
                             splinetime = 0,
                             polytime=3,
                             priors = NULL,
                             intercs = TRUE, 
                             empri = 0.001*nrow(senior.data),
                             ts = "year", 
                             cs = "cowcode", 
                             lags = c("lwelfg", "lln_milex_pc", "lln_gdp", "lln_rents", "lmilex_exp"), 
                             leads = c("lwelfg", "lln_milex_pc", "lln_gdp", "lln_rents", "lmilex_exp"), 
                             logs = NULL, 
                             ords = NULL, 
                             noms = c("juniorcoup_Holger", "seniorcoup_Holger"), 
                             bounds = NULL)

## Run regressions 
z.out.jun <- zelig(juniorcoup_Holger ~ lwelfg + lln_milex_pc + lmilex_exp + leffective_number + 
                   lnopcoups + lln_gdp + lgrowth + lucdp_internal + lethfrac3 + ldpolity2 + lln_rents + lln_instability2 + 
                   junior_spell + junior_spell2 + junior_spell3,
                   data = junior.imputations, model = "probit", robust=T, cluster="cowcode")
summary(z.out.jun)

z.out.sen <- zelig(seniorcoup_Holger ~ lwelfg + lln_milex_pc + lmilex_exp + leffective_number +  
                     lnopcoups + lln_gdp + lgrowth + lucdp_internal + lethfrac3 + ldpolity2 + lln_rents + lln_instability2 + 
                     senior_spell + senior_spell2 + senior_spell3,
                   data = senior.imputations, model = "probit", robust=T, cluster="cowcode")
summary(z.out.sen)
